1   /* Copyright 2002-2016 CS Systèmes d'Information
2    * Licensed to CS Systèmes d'Information (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.estimation.leastsquares;
18  
19  import java.io.BufferedReader;
20  import java.io.File;
21  import java.io.FileInputStream;
22  import java.io.IOException;
23  import java.io.InputStreamReader;
24  import java.io.UnsupportedEncodingException;
25  import java.net.URISyntaxException;
26  import java.text.ParseException;
27  import java.util.ArrayList;
28  import java.util.Collections;
29  import java.util.Comparator;
30  import java.util.HashMap;
31  import java.util.List;
32  import java.util.Map;
33  import java.util.NoSuchElementException;
34  import java.util.SortedSet;
35  import java.util.TreeSet;
36  
37  import org.hipparchus.exception.LocalizedCoreFormats;
38  import org.hipparchus.geometry.euclidean.threed.Vector3D;
39  import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer;
40  import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer;
41  import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
42  import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.Decomposition;
43  import org.hipparchus.stat.descriptive.StreamingStatistics;
44  import org.hipparchus.util.FastMath;
45  import org.hipparchus.util.Precision;
46  import org.junit.Assert;
47  import org.junit.Test;
48  import org.orekit.KeyValueFileParser;
49  import org.orekit.Utils;
50  import org.orekit.bodies.CelestialBody;
51  import org.orekit.bodies.CelestialBodyFactory;
52  import org.orekit.bodies.GeodeticPoint;
53  import org.orekit.bodies.OneAxisEllipsoid;
54  import org.orekit.data.DataProvidersManager;
55  import org.orekit.errors.OrekitException;
56  import org.orekit.errors.OrekitMessages;
57  import org.orekit.estimation.measurements.Angular;
58  import org.orekit.estimation.measurements.Bias;
59  import org.orekit.estimation.measurements.EstimatedMeasurement;
60  import org.orekit.estimation.measurements.GroundStation;
61  import org.orekit.estimation.measurements.ObservedMeasurement;
62  import org.orekit.estimation.measurements.OutlierFilter;
63  import org.orekit.estimation.measurements.PV;
64  import org.orekit.estimation.measurements.Range;
65  import org.orekit.estimation.measurements.RangeRate;
66  import org.orekit.estimation.measurements.modifiers.AngularRadioRefractionModifier;
67  import org.orekit.estimation.measurements.modifiers.RangeTroposphericDelayModifier;
68  import org.orekit.forces.drag.Atmosphere;
69  import org.orekit.forces.drag.DTM2000;
70  import org.orekit.forces.drag.DragForce;
71  import org.orekit.forces.drag.DragSensitive;
72  import org.orekit.forces.drag.IsotropicDrag;
73  import org.orekit.forces.drag.MarshallSolarActivityFutureEstimation;
74  import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
75  import org.orekit.forces.gravity.OceanTides;
76  import org.orekit.forces.gravity.Relativity;
77  import org.orekit.forces.gravity.SolidTides;
78  import org.orekit.forces.gravity.ThirdBodyAttraction;
79  import org.orekit.forces.gravity.potential.GravityFieldFactory;
80  import org.orekit.forces.gravity.potential.ICGEMFormatReader;
81  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
82  import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
83  import org.orekit.forces.radiation.RadiationSensitive;
84  import org.orekit.forces.radiation.SolarRadiationPressure;
85  import org.orekit.frames.Frame;
86  import org.orekit.frames.FramesFactory;
87  import org.orekit.frames.TopocentricFrame;
88  import org.orekit.frames.Transform;
89  import org.orekit.models.AtmosphericRefractionModel;
90  import org.orekit.models.earth.EarthITU453AtmosphereRefraction;
91  import org.orekit.models.earth.SaastamoinenModel;
92  import org.orekit.orbits.CartesianOrbit;
93  import org.orekit.orbits.CircularOrbit;
94  import org.orekit.orbits.EquinoctialOrbit;
95  import org.orekit.orbits.KeplerianOrbit;
96  import org.orekit.orbits.Orbit;
97  import org.orekit.orbits.PositionAngle;
98  import org.orekit.propagation.SpacecraftState;
99  import org.orekit.propagation.analytical.tle.TLE;
100 import org.orekit.propagation.analytical.tle.TLEPropagator;
101 import org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder;
102 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
103 import org.orekit.time.AbsoluteDate;
104 import org.orekit.time.ChronologicalComparator;
105 import org.orekit.time.TimeScalesFactory;
106 import org.orekit.utils.Constants;
107 import org.orekit.utils.IERSConventions;
108 import org.orekit.utils.PVCoordinates;
109 import org.orekit.utils.ParameterDriver;
110 import org.orekit.utils.ParameterDriversList;
111 import org.orekit.utils.ParameterDriversList.DelegatingDriver;
112 
113 public class OrbitDeterminationTest {
114 
115     @Test
116     // Orbit determination for Lageos2 based on SLR (range) measurements
117     public void testLageos2()
118         throws URISyntaxException, IllegalArgumentException, IOException,
119                OrekitException, ParseException {
120 
121         // input in tutorial resources directory/output
122         final String inputPath = OrbitDeterminationTest.class.getClassLoader().getResource("orbit-determination/Lageos2/od_test_Lageos2.in").toURI().getPath();
123         final File input  = new File(inputPath);
124 
125         // configure Orekit data acces
126         Utils.setDataRoot("orbit-determination/Lageos2:potential/icgem-format");
127         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
128 
129         //orbit determination run.
130         ResultOD odLageos2 = run(input);
131 
132         //test
133         //definition of the accuracy for the test
134         final double distanceAccuracy = 0.1;
135         final double velocityAccuracy = 1e-4;
136 
137         //test on the convergence
138         final int numberOfIte  = 4;
139         final int numberOfEval = 4;
140 
141         Assert.assertEquals(numberOfIte, odLageos2.getNumberOfIteration());
142         Assert.assertEquals(numberOfEval, odLageos2.getNumberOfEvaluation());
143 
144         //test on the estimated position and velocity
145         final Vector3D estimatedPos = odLageos2.getEstimatedPV().getPosition();
146         final Vector3D estimatedVel = odLageos2.getEstimatedPV().getVelocity();
147         final Vector3D refPos = new Vector3D(-5532124.989973327, 10025700.01763335, -3578940.840115321);
148         final Vector3D refVel = new Vector3D(-3871.2736402553, -607.8775965705, 4280.9744110925);
149         Assert.assertEquals(0.0, Vector3D.distance(refPos, estimatedPos), distanceAccuracy);
150         Assert.assertEquals(0.0, Vector3D.distance(refVel, estimatedVel), velocityAccuracy);
151 
152         //test on measurements parameters
153         final List<DelegatingDriver> list = new ArrayList<DelegatingDriver>();
154         list.addAll(odLageos2.measurementsParameters.getDrivers());
155         sortParametersChanges(list);
156         final double[] stationOffSet = { -1.351682,  -2.180542,  -5.278784 };
157         final double rangeBias = -7.923393;
158         Assert.assertEquals(stationOffSet[0], list.get(0).getValue(), distanceAccuracy);
159         Assert.assertEquals(stationOffSet[1], list.get(1).getValue(), distanceAccuracy);
160         Assert.assertEquals(stationOffSet[2], list.get(2).getValue(), distanceAccuracy);
161         Assert.assertEquals(rangeBias,        list.get(3).getValue(), distanceAccuracy);
162 
163         //test on statistic for the range residuals
164         final long nbRange = 258;
165         final double[] RefStatRange = { -2.795816, 6.171529, 0.310848, 1.657809 };
166         Assert.assertEquals(nbRange, odLageos2.getRangeStat().getN());
167         Assert.assertEquals(RefStatRange[0], odLageos2.getRangeStat().getMin(),               distanceAccuracy);
168         Assert.assertEquals(RefStatRange[1], odLageos2.getRangeStat().getMax(),               distanceAccuracy);
169         Assert.assertEquals(RefStatRange[2], odLageos2.getRangeStat().getMean(),              distanceAccuracy);
170         Assert.assertEquals(RefStatRange[3], odLageos2.getRangeStat().getStandardDeviation(), distanceAccuracy);
171 
172     }
173 
174     @Test
175     // Orbit determination for range, azimuth elevation measurements
176     public void testW3B()
177         throws URISyntaxException, IllegalArgumentException, IOException,
178               OrekitException, ParseException {
179 
180         // input in tutorial resources directory/output
181         final String inputPath = OrbitDeterminationTest.class.getClassLoader().getResource("orbit-determination/W3B/od_test_W3.in").toURI().getPath();
182         final File input  = new File(inputPath);
183 
184         // configure Orekit data access
185         Utils.setDataRoot("orbit-determination/W3B:potential/icgem-format");
186         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
187 
188         //orbit determination run.
189         ResultOD odsatW3 = run(input);
190 
191         //test
192         //definition of the accuracy for the test
193         final double distanceAccuracy = 1e-1;
194         final double velocityAccuracy = 1e-4;
195         // angle unit is radian
196         final double angleAccuracy = 1e-5;
197         final double dimensionLessCoef = 1e-3;
198 
199         //test on the convergence
200         final int numberOfIte  = 4;
201         final int numberOfEval = 7;
202 
203         Assert.assertEquals(numberOfIte, odsatW3.getNumberOfIteration());
204         Assert.assertEquals(numberOfEval, odsatW3.getNumberOfEvaluation());
205 
206         //test on the estimated position and velocity
207         final Vector3D estimatedPos = odsatW3.getEstimatedPV().getPosition();
208         final Vector3D estimatedVel = odsatW3.getEstimatedPV().getVelocity();
209         final Vector3D refPos = new Vector3D( -40541691.15225173, -9903912.495473776, 208037.5511875451);
210         final Vector3D refVel = new Vector3D( 759.0248098953, -1476.58298279, 54.7065550778);
211         Assert.assertEquals(0.0, Vector3D.distance(refPos, estimatedPos), distanceAccuracy);
212         Assert.assertEquals(0.0, Vector3D.distance(refVel, estimatedVel), velocityAccuracy);
213 
214 
215         //test on propagator parameters
216         final double dragCoef  = 0.215421;
217         final double SRPCoef = 112.336693;
218         Assert.assertEquals(dragCoef, odsatW3.propagatorParameters.getDrivers().get(0).getValue(), 3e-3);
219         Assert.assertEquals(SRPCoef,  odsatW3.propagatorParameters.getDrivers().get(1).getValue(), dimensionLessCoef);
220 
221         //test on measurements parameters
222         final List<DelegatingDriver> list = new ArrayList<DelegatingDriver>();
223         list.addAll(odsatW3.measurementsParameters.getDrivers());
224         sortParametersChanges(list);
225 
226         //station CastleRock
227         final double[] CastleAzElBias = { 0.061500,  -0.004955};
228         final double CastleRangeBias = 11461.679992;
229         Assert.assertEquals(CastleAzElBias[0], FastMath.toDegrees(list.get(0).getValue()), angleAccuracy);
230         Assert.assertEquals(CastleAzElBias[1], FastMath.toDegrees(list.get(1).getValue()), angleAccuracy);
231         Assert.assertEquals(CastleRangeBias,   list.get(2).getValue(),                     distanceAccuracy);
232 
233         //station Fucino
234         final double[] FucAzElBias = {-0.055555,  0.075471};
235         final double FucRangeBias = 13461.362346;
236         Assert.assertEquals(FucAzElBias[0], FastMath.toDegrees(list.get(3).getValue()), angleAccuracy);
237         Assert.assertEquals(FucAzElBias[1], FastMath.toDegrees(list.get(4).getValue()), angleAccuracy);
238         Assert.assertEquals(FucRangeBias,   list.get(5).getValue(),                     distanceAccuracy);
239 
240         //station Kumsan
241         final double[] KumAzElBias = { -0.025227,  -0.054883};
242         final double KumRangeBias = 13521.236722;
243         Assert.assertEquals(KumAzElBias[0], FastMath.toDegrees(list.get(6).getValue()), angleAccuracy);
244         Assert.assertEquals(KumAzElBias[1], FastMath.toDegrees(list.get(7).getValue()), angleAccuracy);
245         Assert.assertEquals(KumRangeBias,   list.get(8).getValue(),                     distanceAccuracy);
246 
247         //station Pretoria
248         final double[] PreAzElBias = { 0.029580,  0.011505};
249         final double PreRangeBias = 13373.739538;
250         Assert.assertEquals(PreAzElBias[0], FastMath.toDegrees(list.get( 9).getValue()), angleAccuracy);
251         Assert.assertEquals(PreAzElBias[1], FastMath.toDegrees(list.get(10).getValue()), angleAccuracy);
252         Assert.assertEquals(PreRangeBias,   list.get(11).getValue(),                     distanceAccuracy);
253 
254         //station Uralla
255         final double[] UraAzElBias = { 0.168236,  -0.121083};
256         final double UraRangeBias = 13294.762426;
257         Assert.assertEquals(UraAzElBias[0], FastMath.toDegrees(list.get(12).getValue()), angleAccuracy);
258         Assert.assertEquals(UraAzElBias[1], FastMath.toDegrees(list.get(13).getValue()), angleAccuracy);
259         Assert.assertEquals(UraRangeBias,   list.get(14).getValue(),                     distanceAccuracy);
260 
261         //test on statistic for the range residuals
262         final long nbRange = 182;
263         //statistics for the range residual (min, max, mean, std)
264         final double[] RefStatRange = { -18.407301888801157, 11.21531879529357, -1.3004374373196126E-5, 4.396892027885725 };
265         Assert.assertEquals(nbRange, odsatW3.getRangeStat().getN());
266         Assert.assertEquals(RefStatRange[0], odsatW3.getRangeStat().getMin(),               distanceAccuracy);
267         Assert.assertEquals(RefStatRange[1], odsatW3.getRangeStat().getMax(),               distanceAccuracy);
268         Assert.assertEquals(RefStatRange[2], odsatW3.getRangeStat().getMean(),              distanceAccuracy);
269         Assert.assertEquals(RefStatRange[3], odsatW3.getRangeStat().getStandardDeviation(), distanceAccuracy);
270 
271         //test on statistic for the azimuth residuals
272         final long nbAzi = 339;
273         //statistics for the azimuth residual (min, max, mean, std)
274         final double[] RefStatAzi = { -0.04334642535664204, 0.025348556908738606, -3.471009794826308E-11, 0.010532646772836074 };
275         Assert.assertEquals(nbAzi, odsatW3.getAzimStat().getN());
276         Assert.assertEquals(RefStatAzi[0], odsatW3.getAzimStat().getMin(),               angleAccuracy);
277         Assert.assertEquals(RefStatAzi[1], odsatW3.getAzimStat().getMax(),               angleAccuracy);
278         Assert.assertEquals(RefStatAzi[2], odsatW3.getAzimStat().getMean(),              angleAccuracy);
279         Assert.assertEquals(RefStatAzi[3], odsatW3.getAzimStat().getStandardDeviation(), angleAccuracy);
280 
281         //test on statistic for the azimuth residuals
282         final long nbEle = 339;
283         //statistics for the azimuth residual (min, max, mean, std)
284         final double[] RefStatEle = { -0.024901037294786422, 0.0549815182281244, -1.2496574162566745E-11, 0.011712542337682996 };
285         Assert.assertEquals(nbEle, odsatW3.getElevStat().getN());
286         Assert.assertEquals(RefStatEle[0], odsatW3.getElevStat().getMin(),               angleAccuracy);
287         Assert.assertEquals(RefStatEle[1], odsatW3.getElevStat().getMax(),               angleAccuracy);
288         Assert.assertEquals(RefStatEle[2], odsatW3.getElevStat().getMean(),              angleAccuracy);
289         Assert.assertEquals(RefStatEle[3], odsatW3.getElevStat().getStandardDeviation(), angleAccuracy);
290 
291     }
292 
293    private class ResultOD {
294        private int numberOfIteration;
295        private int numberOfEvaluation;
296        private PVCoordinates estimatedPV;
297        private StreamingStatistics rangeStat;
298        private StreamingStatistics azimStat;
299        private StreamingStatistics elevStat;
300        private  ParameterDriversList propagatorParameters  ;
301        private  ParameterDriversList measurementsParameters;
302        ResultOD(ParameterDriversList  propagatorParameters,
303                 ParameterDriversList  measurementsParameters,
304                 int numberOfIteration, int numberOfEvaluation, PVCoordinates estimatedPV,
305                 StreamingStatistics rangeStat, StreamingStatistics rangeRateStat,
306                 StreamingStatistics azimStat, StreamingStatistics elevStat,
307                 StreamingStatistics posStat, StreamingStatistics velStat) {
308 
309            this.propagatorParameters   = propagatorParameters;
310            this.measurementsParameters = measurementsParameters;
311            this.numberOfIteration      = numberOfIteration;
312            this.numberOfEvaluation     = numberOfEvaluation;
313            this.estimatedPV            = estimatedPV;
314            this.rangeStat              =  rangeStat;
315            this.azimStat               = azimStat;
316            this.elevStat               = elevStat;
317        }
318 
319 
320     public int getNumberOfIteration() {
321         return numberOfIteration;
322     }
323 
324 
325     public int getNumberOfEvaluation() {
326         return numberOfEvaluation;
327     }
328 
329 
330     public PVCoordinates getEstimatedPV() {
331         return estimatedPV;
332     }
333 
334 
335     public StreamingStatistics getRangeStat() {
336         return rangeStat;
337     }
338 
339 
340 
341     public StreamingStatistics getAzimStat() {
342         return azimStat;
343     }
344 
345 
346     public StreamingStatistics getElevStat() {
347         return elevStat;
348     }
349 
350 
351    }
352 
353     private ResultOD run(final File input)
354         throws IOException, IllegalArgumentException, OrekitException, ParseException {
355 
356         // read input parameters
357         KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
358         parser.parseInput(input.getAbsolutePath(), new FileInputStream(input));
359 
360         // log file
361 
362         final RangeLog     rangeLog     = new RangeLog();
363         final RangeRateLog rangeRateLog = new RangeRateLog();
364         final AzimuthLog   azimuthLog   = new AzimuthLog();
365         final ElevationLog elevationLog = new ElevationLog();
366         final PositionLog  positionLog  = new PositionLog();
367         final VelocityLog  velocityLog  = new VelocityLog();
368 
369         // gravity field
370         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-5c.gfc", true));
371         final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
372 
373 
374         // Orbit initial guess
375         final Orbit initialGuess = createOrbit(parser, gravityField.getMu());
376 
377         // IERS conventions
378         final IERSConventions conventions;
379         if (!parser.containsKey(ParameterKey.IERS_CONVENTIONS)) {
380             conventions = IERSConventions.IERS_2010;
381         } else {
382             conventions = IERSConventions.valueOf("IERS_" + parser.getInt(ParameterKey.IERS_CONVENTIONS));
383         }
384 
385         // central body
386         final OneAxisEllipsoid body = createBody(parser);
387 
388         // propagator builder
389         final NumericalPropagatorBuilder propagatorBuilder =
390                         createPropagatorBuilder(parser, conventions, gravityField, body, initialGuess);
391 
392         // estimator
393         final BatchLSEstimator estimator = createEstimator(parser, propagatorBuilder);
394 
395         // measurements
396         final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
397         for (final String fileName : parser.getStringsList(ParameterKey.MEASUREMENTS_FILES, ',')) {
398             measurements.addAll(readMeasurements(new File(input.getParentFile(), fileName),
399                                                  createStationsData(parser, body),
400                                                  createPVData(parser),
401                                                  createSatRangeBias(parser),
402                                                  createWeights(parser),
403                                                  createRangeOutliersManager(parser),
404                                                  createRangeRateOutliersManager(parser),
405                                                  createAzElOutliersManager(parser),
406                                                  createPVOutliersManager(parser)));
407         }
408         for (ObservedMeasurement<?> measurement : measurements) {
409             estimator.addMeasurement(measurement);
410         }
411 
412         Orbit estimated = estimator.estimate().getInitialState().getOrbit();
413 
414         // compute some statistics
415         for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
416             if (entry.getKey() instanceof Range) {
417                 @SuppressWarnings("unchecked")
418                 EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
419                 rangeLog.add(evaluation);
420             } else if (entry.getKey() instanceof RangeRate) {
421                 @SuppressWarnings("unchecked")
422                 EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
423                 rangeRateLog.add(evaluation);
424             } else if (entry.getKey() instanceof Angular) {
425                 @SuppressWarnings("unchecked")
426                 EstimatedMeasurement<Angular> evaluation = (EstimatedMeasurement<Angular>) entry.getValue();
427                 azimuthLog.add(evaluation);
428                 elevationLog.add(evaluation);
429             } else if (entry.getKey() instanceof PV) {
430                 @SuppressWarnings("unchecked")
431                 EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
432                 positionLog.add(evaluation);
433                 velocityLog.add(evaluation);
434             }
435         }
436 
437         // parmaters and measurements.
438         final ParameterDriversList propagatorParameters   = estimator.getPropagatorParametersDrivers(true);
439         final ParameterDriversList measurementsParameters = estimator.getMeasurementsParametersDrivers(true);
440 
441         //instation of results
442         return new ResultOD(propagatorParameters, measurementsParameters,
443                             estimator.getIterationsCount(), estimator.getEvaluationsCount(), estimated.getPVCoordinates(),
444                             rangeLog.createStatisticsSummary(),  rangeRateLog.createStatisticsSummary(),
445                             azimuthLog.createStatisticsSummary(),  elevationLog.createStatisticsSummary(),
446                             positionLog.createStatisticsSummary(),  velocityLog.createStatisticsSummary());
447     }
448 
449     /** Sort parameters changes.
450      * @param parameters parameters list
451      */
452     private void sortParametersChanges(List<? extends ParameterDriver> parameters) {
453 
454         // sort the parameters lexicographically
455         Collections.sort(parameters, new Comparator<ParameterDriver>() {
456             /** {@inheritDoc} */
457             @Override
458             public int compare(final ParameterDriver pd1, final ParameterDriver pd2) {
459                 return pd1.getName().compareTo(pd2.getName());
460             }
461 
462         });
463     }
464 
465 
466     /** Create a propagator builder from input parameters
467      * @param parser input file parser
468      * @param conventions IERS conventions to use
469      * @param gravityField gravity field
470      * @param body central body
471      * @param orbit first orbit estimate
472      * @return propagator builder
473      * @throws NoSuchElementException if input parameters are missing
474      * @throws OrekitException if body frame cannot be created
475      */
476     private NumericalPropagatorBuilder createPropagatorBuilder(final KeyValueFileParser<ParameterKey> parser,
477                                                                final IERSConventions conventions,
478                                                                final NormalizedSphericalHarmonicsProvider gravityField,
479                                                                final OneAxisEllipsoid body,
480                                                                final Orbit orbit)
481         throws NoSuchElementException, OrekitException {
482 
483         final double minStep;
484         if (!parser.containsKey(ParameterKey.PROPAGATOR_MIN_STEP)) {
485             minStep = 0.001;
486         } else {
487             minStep = parser.getDouble(ParameterKey.PROPAGATOR_MIN_STEP);
488         }
489 
490         final double maxStep;
491         if (!parser.containsKey(ParameterKey.PROPAGATOR_MAX_STEP)) {
492             maxStep = 300;
493         } else {
494             maxStep = parser.getDouble(ParameterKey.PROPAGATOR_MAX_STEP);
495         }
496 
497         final double dP;
498         if (!parser.containsKey(ParameterKey.PROPAGATOR_POSITION_ERROR)) {
499             dP = 10.0;
500         } else {
501             dP = parser.getDouble(ParameterKey.PROPAGATOR_POSITION_ERROR);
502         }
503 
504         final double positionScale;
505         if (!parser.containsKey(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE)) {
506             positionScale = dP;
507         } else {
508             positionScale = parser.getDouble(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE);
509         }
510         final NumericalPropagatorBuilder propagatorBuilder =
511                         new NumericalPropagatorBuilder(orbit,
512                                                        new DormandPrince853IntegratorBuilder(minStep, maxStep, dP),
513                                                        PositionAngle.MEAN,
514                                                        positionScale);
515 
516         // initial mass
517         final double mass;
518         if (!parser.containsKey(ParameterKey.MASS)) {
519             mass = 1000.0;
520         } else {
521             mass = parser.getDouble(ParameterKey.MASS);
522         }
523         propagatorBuilder.setMass(mass);
524 
525         // gravity field force model
526         propagatorBuilder.addForceModel(new HolmesFeatherstoneAttractionModel(body.getBodyFrame(), gravityField));
527 
528         // ocean tides force model
529         if (parser.containsKey(ParameterKey.OCEAN_TIDES_DEGREE) &&
530             parser.containsKey(ParameterKey.OCEAN_TIDES_ORDER)) {
531             final int degree = parser.getInt(ParameterKey.OCEAN_TIDES_DEGREE);
532             final int order  = parser.getInt(ParameterKey.OCEAN_TIDES_ORDER);
533             if (degree > 0 && order > 0) {
534                 propagatorBuilder.addForceModel(new OceanTides(body.getBodyFrame(),
535                                                                gravityField.getAe(), gravityField.getMu(),
536                                                                degree, order, conventions,
537                                                                TimeScalesFactory.getUT1(conventions, true)));
538             }
539         }
540 
541         // solid tides force model
542         List<CelestialBody> solidTidesBodies = new ArrayList<CelestialBody>();
543         if (parser.containsKey(ParameterKey.SOLID_TIDES_SUN) &&
544             parser.getBoolean(ParameterKey.SOLID_TIDES_SUN)) {
545             solidTidesBodies.add(CelestialBodyFactory.getSun());
546         }
547         if (parser.containsKey(ParameterKey.SOLID_TIDES_MOON) &&
548             parser.getBoolean(ParameterKey.SOLID_TIDES_MOON)) {
549             solidTidesBodies.add(CelestialBodyFactory.getMoon());
550         }
551         if (!solidTidesBodies.isEmpty()) {
552             propagatorBuilder.addForceModel(new SolidTides(body.getBodyFrame(),
553                                                            gravityField.getAe(), gravityField.getMu(),
554                                                            gravityField.getTideSystem(), conventions,
555                                                            TimeScalesFactory.getUT1(conventions, true),
556                                                            solidTidesBodies.toArray(new CelestialBody[solidTidesBodies.size()])));
557         }
558 
559         // third body attraction
560         if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) &&
561             parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) {
562             propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
563         }
564         if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) &&
565             parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) {
566             propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
567         }
568 
569         // drag
570         if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) {
571             final double  cd          = parser.getDouble(ParameterKey.DRAG_CD);
572             final double  area        = parser.getDouble(ParameterKey.DRAG_AREA);
573             final boolean cdEstimated = parser.getBoolean(ParameterKey.DRAG_CD_ESTIMATED);
574 
575             MarshallSolarActivityFutureEstimation msafe =
576                             new MarshallSolarActivityFutureEstimation("(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)",
577             MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
578             DataProvidersManager manager = DataProvidersManager.getInstance();
579             manager.feed(msafe.getSupportedNames(), msafe);
580             Atmosphere atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), body);
581             propagatorBuilder.addForceModel(new DragForce(atmosphere, new IsotropicDrag(area, cd)));
582             if (cdEstimated) {
583                 for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
584                     if (driver.getName().equals(DragSensitive.DRAG_COEFFICIENT)) {
585                         driver.setSelected(true);
586                     }
587                 }
588             }
589         }
590 
591         // solar radiation pressure
592         if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) {
593             final double  cr          = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR);
594             final double  area        = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_AREA);
595             final boolean cREstimated = parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE_CR_ESTIMATED);
596 
597             propagatorBuilder.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(),
598                                                                        body.getEquatorialRadius(),
599                                                                        new IsotropicRadiationSingleCoefficient(area, cr)));
600             if (cREstimated) {
601                 for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
602                     if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
603                         driver.setSelected(true);
604                     }
605                 }
606             }
607         }
608 
609         // post-Newtonian correction force due to general relativity
610         if (parser.containsKey(ParameterKey.GENERAL_RELATIVITY) && parser.getBoolean(ParameterKey.GENERAL_RELATIVITY)) {
611             propagatorBuilder.addForceModel(new Relativity(gravityField.getMu()));
612         }
613 
614         return propagatorBuilder;
615 
616     }
617 
618     /** Create a gravity field from input parameters
619      * @param parser input file parser
620      * @return gravity field
621      * @throws NoSuchElementException if input parameters are missing
622      * @throws OrekitException if body frame cannot be created
623      */
624     private NormalizedSphericalHarmonicsProvider createGravityField(final KeyValueFileParser<ParameterKey> parser)
625         throws NoSuchElementException, OrekitException {
626         final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE);
627         final int order  = FastMath.min(degree, parser.getInt(ParameterKey.CENTRAL_BODY_ORDER));
628         return GravityFieldFactory.getNormalizedProvider(degree, order);
629     }
630 
631     /** Create an orbit from input parameters
632      * @param parser input file parser
633      * @param mu     central attraction coefficient
634      * @throws NoSuchElementException if input parameters are missing
635      * @throws OrekitException if body frame cannot be created
636      */
637     private OneAxisEllipsoid createBody(final KeyValueFileParser<ParameterKey> parser)
638         throws NoSuchElementException, OrekitException {
639 
640         final Frame bodyFrame;
641         if (!parser.containsKey(ParameterKey.BODY_FRAME)) {
642             bodyFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
643         } else {
644             bodyFrame = parser.getEarthFrame(ParameterKey.BODY_FRAME);
645         }
646 
647         final double equatorialRadius;
648         if (!parser.containsKey(ParameterKey.BODY_EQUATORIAL_RADIUS)) {
649             equatorialRadius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
650         } else {
651             equatorialRadius = parser.getDouble(ParameterKey.BODY_EQUATORIAL_RADIUS);
652         }
653 
654         final double flattening;
655         if (!parser.containsKey(ParameterKey.BODY_INVERSE_FLATTENING)) {
656             flattening = Constants.WGS84_EARTH_FLATTENING;
657         } else {
658             flattening = 1.0 / parser.getDouble(ParameterKey.BODY_INVERSE_FLATTENING);
659         }
660 
661         return new OneAxisEllipsoid(equatorialRadius, flattening, bodyFrame);
662 
663     }
664 
665     /** Create an orbit from input parameters
666      * @param parser input file parser
667      * @param mu     central attraction coefficient
668      * @throws NoSuchElementException if input parameters are missing
669      * @throws OrekitException if inertial frame cannot be created
670      */
671     private Orbit createOrbit(final KeyValueFileParser<ParameterKey> parser,
672                               final double mu)
673         throws NoSuchElementException, OrekitException {
674 
675         final Frame frame;
676         if (!parser.containsKey(ParameterKey.INERTIAL_FRAME)) {
677             frame = FramesFactory.getEME2000();
678         } else {
679             frame = parser.getInertialFrame(ParameterKey.INERTIAL_FRAME);
680         }
681 
682         // Orbit definition
683         PositionAngle angleType = PositionAngle.MEAN;
684         if (parser.containsKey(ParameterKey.ORBIT_ANGLE_TYPE)) {
685             angleType = PositionAngle.valueOf(parser.getString(ParameterKey.ORBIT_ANGLE_TYPE).toUpperCase());
686         }
687         if (parser.containsKey(ParameterKey.ORBIT_KEPLERIAN_A)) {
688             return new KeplerianOrbit(parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_A),
689                                       parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_E),
690                                       parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_I),
691                                       parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_PA),
692                                       parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_RAAN),
693                                       parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_ANOMALY),
694                                       angleType,
695                                       frame,
696                                       parser.getDate(ParameterKey.ORBIT_DATE,
697                                                      TimeScalesFactory.getUTC()),
698                                       mu);
699         } else if (parser.containsKey(ParameterKey.ORBIT_EQUINOCTIAL_A)) {
700             return new EquinoctialOrbit(parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_A),
701                                         parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EX),
702                                         parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EY),
703                                         parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HX),
704                                         parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HY),
705                                         parser.getAngle(ParameterKey.ORBIT_EQUINOCTIAL_LAMBDA),
706                                         angleType,
707                                         frame,
708                                         parser.getDate(ParameterKey.ORBIT_DATE,
709                                                        TimeScalesFactory.getUTC()),
710                                         mu);
711         } else if (parser.containsKey(ParameterKey.ORBIT_CIRCULAR_A)) {
712             return new CircularOrbit(parser.getDouble(ParameterKey.ORBIT_CIRCULAR_A),
713                                      parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EX),
714                                      parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EY),
715                                      parser.getAngle(ParameterKey.ORBIT_CIRCULAR_I),
716                                      parser.getAngle(ParameterKey.ORBIT_CIRCULAR_RAAN),
717                                      parser.getAngle(ParameterKey.ORBIT_CIRCULAR_ALPHA),
718                                      angleType,
719                                      frame,
720                                      parser.getDate(ParameterKey.ORBIT_DATE,
721                                                     TimeScalesFactory.getUTC()),
722                                      mu);
723         } else if (parser.containsKey(ParameterKey.ORBIT_TLE_LINE_1)) {
724             final String line1 = parser.getString(ParameterKey.ORBIT_TLE_LINE_1);
725             final String line2 = parser.getString(ParameterKey.ORBIT_TLE_LINE_2);
726             final TLE tle = new TLE(line1, line2);
727 
728             TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
729            // propagator.setEphemerisMode();
730 
731             AbsoluteDate initDate = tle.getDate();
732             SpacecraftState initialState = propagator.getInitialState();
733 
734 
735             //Transformation from TEME to frame.
736             Transform t =FramesFactory.getTEME().getTransformTo(FramesFactory.getEME2000(), initDate.getDate());
737             return new CartesianOrbit( t.transformPVCoordinates(initialState.getPVCoordinates()) ,
738                                       frame,
739                                       initDate,
740                                        mu);
741 
742 
743         } else {
744             final double[] pos = {parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PX),
745                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PY),
746                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PZ)};
747             final double[] vel = {parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VX),
748                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VY),
749                                   parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VZ)};
750 
751             return new CartesianOrbit(new PVCoordinates(new Vector3D(pos), new Vector3D(vel)),
752                                       frame,
753                                       parser.getDate(ParameterKey.ORBIT_DATE,
754                                                      TimeScalesFactory.getUTC()),
755                                       mu);
756         }
757 
758     }
759 
760     /** Set up range bias due to transponder delay.
761      * @param parser input file parser
762      * @param range bias (may be null if bias is fixed to zero)
763      * @exception OrekitException if bias initial value cannot be set
764      */
765     private Bias<Range> createSatRangeBias(final KeyValueFileParser<ParameterKey> parser)
766         throws OrekitException {
767 
768         // transponder delay
769         final double transponderDelayBias;
770         if (!parser.containsKey(ParameterKey.TRANSPONDER_DELAY_BIAS)) {
771             transponderDelayBias = 0;
772         } else {
773             transponderDelayBias = parser.getDouble(ParameterKey.TRANSPONDER_DELAY_BIAS);
774         }
775 
776         final double transponderDelayBiasMin;
777         if (!parser.containsKey(ParameterKey.TRANSPONDER_DELAY_BIAS_MIN)) {
778             transponderDelayBiasMin = Double.NEGATIVE_INFINITY;
779         } else {
780             transponderDelayBiasMin = parser.getDouble(ParameterKey.TRANSPONDER_DELAY_BIAS_MIN);
781         }
782 
783         final double transponderDelayBiasMax;
784         if (!parser.containsKey(ParameterKey.TRANSPONDER_DELAY_BIAS_MAX)) {
785             transponderDelayBiasMax = Double.NEGATIVE_INFINITY;
786         } else {
787             transponderDelayBiasMax = parser.getDouble(ParameterKey.TRANSPONDER_DELAY_BIAS_MAX);
788         }
789 
790         // bias estimation flag
791         final boolean transponderDelayBiasEstimated;
792         if (!parser.containsKey(ParameterKey.TRANSPONDER_DELAY_BIAS_ESTIMATED)) {
793             transponderDelayBiasEstimated = false;
794         } else {
795             transponderDelayBiasEstimated = parser.getBoolean(ParameterKey.TRANSPONDER_DELAY_BIAS_ESTIMATED);
796         }
797 
798         if (FastMath.abs(transponderDelayBias) >= Precision.SAFE_MIN || transponderDelayBiasEstimated) {
799             // bias is either non-zero or will be estimated,
800             // we really need to create a modifier for this
801             final Bias<Range> bias = new Bias<Range>(new String[] {
802                                                          "transponder delay bias",
803                                                      },
804                                                      new double[] {
805                                                          transponderDelayBias             
806                                                      },
807                                                      new double[] {
808                                                          1.0             
809                                                      },
810                                                      new double[] {
811                                                          transponderDelayBiasMin
812                                                      },
813                                                      new double[] {
814                                                          transponderDelayBiasMax
815                                                      });
816             bias.getParametersDrivers().get(0).setSelected(transponderDelayBiasEstimated);
817             return bias;
818         } else {
819             // fixed zero bias, we don't need any modifier
820             return null;
821         }
822 
823     }
824 
825     /** Set up stations.
826      * @param parser input file parser
827      * @param body central body
828      * @return name to station data map
829      * @exception OrekitException if some frame transforms cannot be computed
830      * @throws NoSuchElementException if input parameters are missing
831      */
832     private Map<String, StationData> createStationsData(final KeyValueFileParser<ParameterKey> parser,
833                                                         final OneAxisEllipsoid body)
834         throws OrekitException, NoSuchElementException {
835 
836         final Map<String, StationData> stations       = new HashMap<String, StationData>();
837 
838         final String[]  stationNames                  = parser.getStringArray(ParameterKey.GROUND_STATION_NAME);
839         final double[]  stationLatitudes              = parser.getAngleArray(ParameterKey.GROUND_STATION_LATITUDE);
840         final double[]  stationLongitudes             = parser.getAngleArray(ParameterKey.GROUND_STATION_LONGITUDE);
841         final double[]  stationAltitudes              = parser.getDoubleArray(ParameterKey.GROUND_STATION_ALTITUDE);
842         final boolean[] stationPositionEstimated      = parser.getBooleanArray(ParameterKey.GROUND_STATION_POSITION_ESTIMATED);
843         final double[]  stationRangeSigma             = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_SIGMA);
844         final double[]  stationRangeBias              = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_BIAS);
845         final double[]  stationRangeBiasMin           = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_BIAS_MIN);
846         final double[]  stationRangeBiasMax           = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_BIAS_MAX);
847         final boolean[] stationRangeBiasEstimated     = parser.getBooleanArray(ParameterKey.GROUND_STATION_RANGE_BIAS_ESTIMATED);
848         final double[]  stationRangeRateSigma         = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_RATE_SIGMA);
849         final double[]  stationRangeRateBias          = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_RATE_BIAS);
850         final double[]  stationRangeRateBiasMin       = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_RATE_BIAS_MIN);
851         final double[]  stationRangeRateBiasMax       = parser.getDoubleArray(ParameterKey.GROUND_STATION_RANGE_RATE_BIAS_MAX);
852         final boolean[] stationRangeRateBiasEstimated = parser.getBooleanArray(ParameterKey.GROUND_STATION_RANGE_RATE_BIAS_ESTIMATED);
853         final double[]  stationAzimuthSigma           = parser.getAngleArray(ParameterKey.GROUND_STATION_AZIMUTH_SIGMA);
854         final double[]  stationAzimuthBias            = parser.getAngleArray(ParameterKey.GROUND_STATION_AZIMUTH_BIAS);
855         final double[]  stationAzimuthBiasMin         = parser.getAngleArray(ParameterKey.GROUND_STATION_AZIMUTH_BIAS_MIN);
856         final double[]  stationAzimuthBiasMax         = parser.getAngleArray(ParameterKey.GROUND_STATION_AZIMUTH_BIAS_MAX);
857         final double[]  stationElevationSigma         = parser.getAngleArray(ParameterKey.GROUND_STATION_ELEVATION_SIGMA);
858         final double[]  stationElevationBias          = parser.getAngleArray(ParameterKey.GROUND_STATION_ELEVATION_BIAS);
859         final double[]  stationElevationBiasMin       = parser.getAngleArray(ParameterKey.GROUND_STATION_ELEVATION_BIAS_MIN);
860         final double[]  stationElevationBiasMax       = parser.getAngleArray(ParameterKey.GROUND_STATION_ELEVATION_BIAS_MAX);
861         final boolean[] stationAzElBiasesEstimated    = parser.getBooleanArray(ParameterKey.GROUND_STATION_AZ_EL_BIASES_ESTIMATED);
862         final boolean[] stationElevationRefraction    = parser.getBooleanArray(ParameterKey.GROUND_STATION_ELEVATION_REFRACTION_CORRECTION);
863         final boolean[] stationRangeTropospheric      = parser.getBooleanArray(ParameterKey.GROUND_STATION_RANGE_TROPOSPHERIC_CORRECTION);
864         //final boolean[] stationIonosphericCorrection    = parser.getBooleanArray(ParameterKey.GROUND_STATION_IONOSPHERIC_CORRECTION);
865 
866         for (int i = 0; i < stationNames.length; ++i) {
867 
868             // the station itself
869             final GeodeticPoint position = new GeodeticPoint(stationLatitudes[i],
870                                                              stationLongitudes[i],
871                                                              stationAltitudes[i]);
872             final TopocentricFrame topo = new TopocentricFrame(body, position, stationNames[i]);
873             final GroundStation station = new GroundStation(topo);
874             station.getEastOffsetDriver().setSelected(stationPositionEstimated[i]);
875             station.getNorthOffsetDriver().setSelected(stationPositionEstimated[i]);
876             station.getZenithOffsetDriver().setSelected(stationPositionEstimated[i]);
877 
878             // range
879             final double rangeSigma = stationRangeSigma[i];
880             final Bias<Range> rangeBias;
881             if (FastMath.abs(stationRangeBias[i])   >= Precision.SAFE_MIN || stationRangeBiasEstimated[i]) {
882                  rangeBias = new Bias<Range>(new String[] {
883                                                  stationNames[i] + "/range bias",
884                                              },
885                                              new double[] {
886                                                  stationRangeBias[i]
887                                              },
888                                              new double[] {
889                                                  rangeSigma
890                                              },
891                                              new double[] {
892                                                  stationRangeBiasMin[i]
893                                              },
894                                              new double[] {
895                                                  stationRangeBiasMax[i]
896                                              });
897                  rangeBias.getParametersDrivers().get(0).setSelected(stationRangeBiasEstimated[i]);
898             } else {
899                 // bias fixed to zero, we don't need to create a modifier for this
900                 rangeBias = null;
901             }
902 
903             // range rate
904             final double rangeRateSigma = stationRangeRateSigma[i];
905             final Bias<RangeRate> rangeRateBias;
906             if (FastMath.abs(stationRangeRateBias[i])   >= Precision.SAFE_MIN || stationRangeRateBiasEstimated[i]) {
907                 rangeRateBias = new Bias<RangeRate>(new String[] {
908                                                         stationNames[i] + "/range rate bias"
909                                                     },
910                                                     new double[] {
911                                                         stationRangeRateBias[i]
912                                                     },
913                                                     new double[] {
914                                                         rangeRateSigma
915                                                     },
916                                                     new double[] {
917                                                         stationRangeRateBiasMin[i]
918                                                     },
919                                                     new double[] {
920                                                         stationRangeRateBiasMax[i]
921                                                     });
922                 rangeRateBias.getParametersDrivers().get(0).setSelected(stationRangeRateBiasEstimated[i]);
923             } else {
924                 // bias fixed to zero, we don't need to create a modifier for this
925                 rangeRateBias = null;
926             }
927 
928             // angular biases
929             final double[] azELSigma = new double[] {
930                 stationAzimuthSigma[i], stationElevationSigma[i]
931             };
932             final Bias<Angular> azELBias;
933             if (FastMath.abs(stationAzimuthBias[i])   >= Precision.SAFE_MIN ||
934                 FastMath.abs(stationElevationBias[i]) >= Precision.SAFE_MIN ||
935                 stationAzElBiasesEstimated[i]) {
936                 azELBias = new Bias<Angular>(new String[] {
937                                                  stationNames[i] + "/az bias",
938                                                  stationNames[i] + "/el bias"
939                                              },
940                                              new double[] {
941                                                  stationAzimuthBias[i],
942                                                  stationElevationBias[i]             
943                                              },
944                                              azELSigma,
945                                              new double[] {
946                                                  stationAzimuthBiasMin[i],
947                                                  stationElevationBiasMin[i]
948                                              },
949                                              new double[] {
950                                                  stationAzimuthBiasMax[i],
951                                                  stationElevationBiasMax[i]
952                                              });
953                 azELBias.getParametersDrivers().get(0).setSelected(stationAzElBiasesEstimated[i]);
954                 azELBias.getParametersDrivers().get(1).setSelected(stationAzElBiasesEstimated[i]);
955             } else {
956                 // bias fixed to zero, we don't need to create a modifier for this
957                 azELBias = null;
958             }
959 
960             //Refraction correction
961             final AngularRadioRefractionModifier refractionCorrection;
962             if (stationElevationRefraction[i]) {
963                 final double                     altitude        = station.getBaseFrame().getPoint().getAltitude();
964                 final AtmosphericRefractionModel refractionModel = new EarthITU453AtmosphereRefraction(altitude);
965                 refractionCorrection = new AngularRadioRefractionModifier(refractionModel);
966             } else {
967                 refractionCorrection = null;
968             }
969 
970 
971             //Tropospheric correction
972             final RangeTroposphericDelayModifier rangeTroposphericCorrection;
973             if (stationRangeTropospheric[i]) {
974 
975                 final SaastamoinenModel troposphericModel = SaastamoinenModel.getStandardModel();
976 
977                 rangeTroposphericCorrection = new  RangeTroposphericDelayModifier(troposphericModel);
978             } else {
979                 rangeTroposphericCorrection = null;
980             }
981 
982 
983         stations.put(stationNames[i], new StationData(station,
984                                                       rangeSigma,     rangeBias,
985                                                       rangeRateSigma, rangeRateBias,
986                                                       azELSigma,      azELBias,
987                                                       refractionCorrection, rangeTroposphericCorrection));
988         }
989         return stations;
990 
991     }
992 
993     /** Set up weights.
994      * @param parser input file parser
995      * @return base weights
996      * @throws NoSuchElementException if input parameters are missing
997      */
998     private Weights createWeights(final KeyValueFileParser<ParameterKey> parser)
999         throws NoSuchElementException {
1000         return new Weights(parser.getDouble(ParameterKey.RANGE_MEASUREMENTS_BASE_WEIGHT),
1001                            parser.getDouble(ParameterKey.RANGE_RATE_MEASUREMENTS_BASE_WEIGHT),
1002                            new double[] {
1003                                parser.getAngle(ParameterKey.AZIMUTH_MEASUREMENTS_BASE_WEIGHT),
1004                                parser.getAngle(ParameterKey.ELEVATION_MEASUREMENTS_BASE_WEIGHT)
1005                            },
1006                            parser.getDouble(ParameterKey.PV_MEASUREMENTS_BASE_WEIGHT));
1007     }
1008 
1009     /** Set up outliers manager for range measurements.
1010      * @param parser input file parser
1011      * @return outliers manager (null if none configured)
1012      * @throws OrekitException if outliers are partly configured
1013      */
1014     private OutlierFilter<Range> createRangeOutliersManager(final KeyValueFileParser<ParameterKey> parser)
1015         throws OrekitException {
1016         if (parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER) !=
1017             parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION)) {
1018             throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1019                                       ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER.toString().toLowerCase().replace('_', '.') +
1020                                       " and  " +
1021                                       ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION.toString().toLowerCase().replace('_', '.') +
1022                                       " must be both present or both absent");
1023         }
1024         return new OutlierFilter<Range>(parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION),
1025                                         parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER));
1026     }
1027 
1028     /** Set up outliers manager for range-rate measurements.
1029      * @param parser input file parser
1030      * @return outliers manager (null if none configured)
1031      * @throws OrekitException if outliers are partly configured
1032      */
1033     private OutlierFilter<RangeRate> createRangeRateOutliersManager(final KeyValueFileParser<ParameterKey> parser)
1034         throws OrekitException {
1035         if (parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER) !=
1036             parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION)) {
1037             throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1038                                       ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER.toString().toLowerCase().replace('_', '.') +
1039                                       " and  " +
1040                                       ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION.toString().toLowerCase().replace('_', '.') +
1041                                       " must be both present or both absent");
1042         }
1043         return new OutlierFilter<RangeRate>(parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION),
1044                                             parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER));
1045     }
1046 
1047     /** Set up outliers manager for azimuth-elevation measurements.
1048      * @param parser input file parser
1049      * @return outliers manager (null if none configured)
1050      * @throws OrekitException if outliers are partly configured
1051      */
1052     private OutlierFilter<Angular> createAzElOutliersManager(final KeyValueFileParser<ParameterKey> parser)
1053         throws OrekitException {
1054         if (parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER) !=
1055             parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION)) {
1056             throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1057                                       ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER.toString().toLowerCase().replace('_', '.') +
1058                                       " and  " +
1059                                       ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION.toString().toLowerCase().replace('_', '.') +
1060                                       " must be both present or both absent");
1061         }
1062         return new OutlierFilter<Angular>(parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION),
1063                                           parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER));
1064     }
1065 
1066     /** Set up outliers manager for PV measurements.
1067      * @param parser input file parser
1068      * @return outliers manager (null if none configured)
1069      * @throws OrekitException if outliers are partly configured
1070      */
1071     private OutlierFilter<PV> createPVOutliersManager(final KeyValueFileParser<ParameterKey> parser)
1072         throws OrekitException {
1073         if (parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER) !=
1074             parser.containsKey(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION)) {
1075             throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1076                                       ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER.toString().toLowerCase().replace('_', '.') +
1077                                       " and  " +
1078                                       ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION.toString().toLowerCase().replace('_', '.') +
1079                                       " must be both present or both absent");
1080         }
1081         return new OutlierFilter<PV>(parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_STARTING_ITERATION),
1082                                      parser.getInt(ParameterKey.RANGE_OUTLIER_REJECTION_MULTIPLIER));
1083     }
1084 
1085     /** Set up PV data.
1086      * @param parser input file parser
1087      * @return PV data
1088      * @throws NoSuchElementException if input parameters are missing
1089      */
1090     private PVData createPVData(final KeyValueFileParser<ParameterKey> parser)
1091         throws NoSuchElementException {
1092         return new PVData(parser.getDouble(ParameterKey.PV_MEASUREMENTS_POSITION_SIGMA),
1093                           parser.getDouble(ParameterKey.PV_MEASUREMENTS_VELOCITY_SIGMA));
1094     }
1095 
1096     /** Set up estimator.
1097      * @param parser input file parser
1098      * @param propagatorBuilder propagator builder
1099      * @return estimator
1100      * @throws NoSuchElementException if input parameters are missing
1101      * @throws OrekitException if some propagator parameters cannot be retrieved
1102      */
1103     private BatchLSEstimator createEstimator(final KeyValueFileParser<ParameterKey> parser,
1104                                              final NumericalPropagatorBuilder propagatorBuilder)
1105         throws NoSuchElementException, OrekitException {
1106         final boolean optimizerIsLevenbergMarquardt;
1107         if (! parser.containsKey(ParameterKey.ESTIMATOR_OPTIMIZATION_ENGINE)) {
1108             optimizerIsLevenbergMarquardt = true;
1109         } else {
1110             final String engine = parser.getString(ParameterKey.ESTIMATOR_OPTIMIZATION_ENGINE);
1111             optimizerIsLevenbergMarquardt = engine.toLowerCase().contains("levenberg");
1112         }
1113         final LeastSquaresOptimizer optimizer;
1114 
1115         if (optimizerIsLevenbergMarquardt) {
1116             // we want to use a Levenberg-Marquardt optimization engine
1117             final double initialStepBoundFactor;
1118             if (! parser.containsKey(ParameterKey.ESTIMATOR_LEVENBERG_MARQUARDT_INITIAL_STEP_BOUND_FACTOR)) {
1119                 initialStepBoundFactor = 100.0;
1120             } else {
1121                 initialStepBoundFactor = parser.getDouble(ParameterKey.ESTIMATOR_LEVENBERG_MARQUARDT_INITIAL_STEP_BOUND_FACTOR);
1122             }
1123 
1124             optimizer = new LevenbergMarquardtOptimizer().withInitialStepBoundFactor(initialStepBoundFactor);
1125         } else {
1126             // we want to use a Gauss-Newton optimization engine
1127             optimizer = new GaussNewtonOptimizer(Decomposition.QR);
1128         }
1129 
1130         final double convergenceThreshold;
1131         if (! parser.containsKey(ParameterKey.ESTIMATOR_NORMALIZED_PARAMETERS_CONVERGENCE_THRESHOLD)) {
1132             convergenceThreshold = 1.0e-3;
1133         } else {
1134             convergenceThreshold = parser.getDouble(ParameterKey.ESTIMATOR_NORMALIZED_PARAMETERS_CONVERGENCE_THRESHOLD);
1135         }
1136         final int maxIterations;
1137         if (! parser.containsKey(ParameterKey.ESTIMATOR_MAX_ITERATIONS)) {
1138             maxIterations = 10;
1139         } else {
1140             maxIterations = parser.getInt(ParameterKey.ESTIMATOR_MAX_ITERATIONS);
1141         }
1142         final int maxEvaluations;
1143         if (! parser.containsKey(ParameterKey.ESTIMATOR_MAX_EVALUATIONS)) {
1144             maxEvaluations = 20;
1145         } else {
1146             maxEvaluations = parser.getInt(ParameterKey.ESTIMATOR_MAX_EVALUATIONS);
1147         }
1148 
1149         final BatchLSEstimator estimator = new BatchLSEstimator(propagatorBuilder, optimizer);
1150         estimator.setParametersConvergenceThreshold(convergenceThreshold);
1151         estimator.setMaxIterations(maxIterations);
1152         estimator.setMaxEvaluations(maxEvaluations);
1153 
1154         return estimator;
1155 
1156     }
1157 
1158     /** Read a measurements file.
1159      * @param file measurements file
1160      * @param stations name to stations data map
1161      * @param pvData PV measurements data
1162      * @param satRangeBias range bias due to transponder delay
1163      * @param weights base weights for measurements
1164      * @param rangeOutliersManager manager for range measurements outliers (null if none configured)
1165      * @param rangeRateOutliersManager manager for range-rate measurements outliers (null if none configured)
1166      * @param azElOutliersManager manager for azimuth-elevation measurements outliers (null if none configured)
1167      * @param pvOutliersManager manager for PV measurements outliers (null if none configured)
1168      * @return measurements list
1169      */
1170     private List<ObservedMeasurement<?>> readMeasurements(final File file,
1171                                                   final Map<String, StationData> stations,
1172                                                   final PVData pvData,
1173                                                   final Bias<Range> satRangeBias,
1174                                                   final Weights weights,
1175                                                   final OutlierFilter<Range> rangeOutliersManager,
1176                                                   final OutlierFilter<RangeRate> rangeRateOutliersManager,
1177                                                   final OutlierFilter<Angular> azElOutliersManager,
1178                                                   final OutlierFilter<PV> pvOutliersManager)
1179         throws UnsupportedEncodingException, IOException, OrekitException {
1180 
1181         final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
1182         BufferedReader br = null;
1183         try {
1184             br = new BufferedReader(new InputStreamReader(new FileInputStream(file), "UTF-8"));
1185             int lineNumber = 0;
1186             for (String line = br.readLine(); line != null; line = br.readLine()) {
1187                 ++lineNumber;
1188                 line = line.trim();
1189                 if (line.length() > 0 && !line.startsWith("#")) {
1190                     String[] fields = line.split("\\s+");
1191                     if (fields.length < 2) {
1192                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
1193                                                   lineNumber, file.getName(), line);
1194                     }
1195                     switch (fields[1]) {
1196                         case "RANGE" :
1197                             final Range range = new RangeParser().parseFields(fields, stations, pvData,
1198                                                                               satRangeBias, weights,
1199                                                                               line, lineNumber, file.getName());
1200                             if (rangeOutliersManager != null) {
1201                                 range.addModifier(rangeOutliersManager);
1202                             }
1203                             addIfNonZeroWeight(range, measurements);
1204                             break;
1205                         case "RANGE_RATE" :
1206                             final RangeRate rangeRate = new RangeRateParser().parseFields(fields, stations, pvData,
1207                                                                                           satRangeBias, weights,
1208                                                                                           line, lineNumber, file.getName());
1209                             if (rangeOutliersManager != null) {
1210                                 rangeRate.addModifier(rangeRateOutliersManager);
1211                             }
1212                             addIfNonZeroWeight(rangeRate, measurements);
1213                             break;
1214                         case "AZ_EL" :
1215                             final Angular angular = new AzElParser().parseFields(fields, stations, pvData,
1216                                                                                  satRangeBias, weights,
1217                                                                                  line, lineNumber, file.getName());
1218                             if (azElOutliersManager != null) {
1219                                 angular.addModifier(azElOutliersManager);
1220                             }
1221                             addIfNonZeroWeight(angular, measurements);
1222                             break;
1223                         case "PV" :
1224                             final PV pv = new PVParser().parseFields(fields, stations, pvData,
1225                                                                      satRangeBias, weights,
1226                                                                      line, lineNumber, file.getName());
1227                             if (pvOutliersManager != null) {
1228                                 pv.addModifier(pvOutliersManager);
1229                             }
1230                             addIfNonZeroWeight(pv, measurements);
1231                             break;
1232                         default :
1233                             throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1234                                                       "unknown measurement type " + fields[1] +
1235                                                       " at line " + lineNumber +
1236                                                       " in file " + file.getName());
1237                     }
1238                 }
1239             }
1240         } finally {
1241             if (br != null) {
1242                 br.close();
1243             }
1244         }
1245 
1246         if (measurements.isEmpty()) {
1247             throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1248                                       "not measurements read from file " + file.getAbsolutePath());
1249         }
1250 
1251         return measurements;
1252 
1253     }
1254 
1255     /** Add a measurement to a list if it has non-zero weight.
1256      * @param measurement measurement to add
1257      * @param measurements measurements list
1258      */
1259     private void addIfNonZeroWeight(final ObservedMeasurement<?> measurement, final List<ObservedMeasurement<?>> measurements) {
1260         double sum = 0;
1261         for (double w : measurement.getBaseWeight()) {
1262             sum += FastMath.abs(w);
1263         }
1264         if (sum > Precision.SAFE_MIN) {
1265             // we only consider measurements with non-zero weight
1266             measurements.add(measurement);
1267         }
1268     }
1269 
1270     /** Container for stations-related data. */
1271     private static class StationData {
1272 
1273         /** Ground station. */
1274         private final GroundStation station;
1275 
1276         /** Range sigma. */
1277         private final double rangeSigma;
1278 
1279         /** Range bias (may be if bias is fixed to zero). */
1280         private final Bias<Range> rangeBias;
1281 
1282         /** Range rate sigma. */
1283         private final double rangeRateSigma;
1284 
1285         /** Range rate bias (may be null if bias is fixed to zero). */
1286         private final Bias<RangeRate> rangeRateBias;
1287 
1288         /** Azimuth-elevation sigma. */
1289         private final double[] azElSigma;
1290 
1291         /** Azimuth-elevation bias (may be null if bias is fixed to zero). */
1292         private final Bias<Angular> azELBias;
1293 
1294         /** Elevation refraction correction (may be null). */
1295         private final AngularRadioRefractionModifier refractionCorrection;
1296 
1297         /** Tropospheric correction (may be null). */
1298         private final RangeTroposphericDelayModifier rangeTroposphericCorrection;
1299 
1300         /** Simple constructor.
1301          * @param station ground station
1302          * @param rangeSigma range sigma
1303          * @param rangeBias range bias (may be null if bias is fixed to zero)
1304          * @param rangeRateSigma range rate sigma
1305          * @param rangeRateBias range rate bias (may be null if bias is fixed to zero)
1306          * @param azElSigma azimuth-elevation sigma
1307          * @param azELBias azimuth-elevation bias (may be null if bias is fixed to zero)
1308          * @param refractionCorrection refraction correction for elevation (may be null)
1309          * @param rangeTroposphericCorrection tropospheric correction  for the range (may be null)
1310          */
1311         public StationData(final GroundStation station,
1312                            final double rangeSigma, final Bias<Range> rangeBias,
1313                            final double rangeRateSigma, final Bias<RangeRate> rangeRateBias,
1314                            final double[] azElSigma, final Bias<Angular> azELBias,
1315                            final AngularRadioRefractionModifier refractionCorrection,
1316                            final RangeTroposphericDelayModifier rangeTroposphericCorrection) {
1317             this.station                     = station;
1318             this.rangeSigma                  = rangeSigma;
1319             this.rangeBias                   = rangeBias;
1320             this.rangeRateSigma              = rangeRateSigma;
1321             this.rangeRateBias               = rangeRateBias;
1322             this.azElSigma                   = azElSigma.clone();
1323             this.azELBias                    = azELBias;
1324             this.refractionCorrection        = refractionCorrection;
1325             this.rangeTroposphericCorrection = rangeTroposphericCorrection;
1326         }
1327 
1328     }
1329 
1330     /** Container for base weights. */
1331     private static class Weights {
1332 
1333         /** Base weight for range measurements. */
1334         private final double rangeBaseWeight;
1335 
1336         /** Base weight for range rate measurements. */
1337         private final double rangeRateBaseWeight;
1338 
1339         /** Base weight for azimuth-elevation measurements. */
1340         private final double[] azElBaseWeight;
1341 
1342         /** Base weight for PV measurements. */
1343         private final double pvBaseWeight;
1344 
1345         /** Simple constructor.
1346          * @param rangeBaseWeight base weight for range measurements
1347          * @param rangeRateBaseWeight base weight for range rate measurements
1348          * @param azElBaseWeight base weight for azimuth-elevation measurements
1349          * @param pvBaseWeight base weight for PV measurements
1350          */
1351         public Weights(final double rangeBaseWeight,
1352                        final double rangeRateBaseWeight,
1353                        final double[] azElBaseWeight,
1354                        final double pvBaseWeight) {
1355             this.rangeBaseWeight     = rangeBaseWeight;
1356             this.rangeRateBaseWeight = rangeRateBaseWeight;
1357             this.azElBaseWeight      = azElBaseWeight.clone();
1358             this.pvBaseWeight        = pvBaseWeight;
1359         }
1360 
1361     }
1362 
1363     /** Container for Position-velocity data. */
1364     private static class PVData {
1365 
1366         /** Position sigma. */
1367         private final double positionSigma;
1368 
1369         /** Velocity sigma. */
1370         private final double velocitySigma;
1371 
1372         /** Simple constructor.
1373          * @param positionSigma position sigma
1374          * @param velocitySigma velocity sigma
1375          */
1376         public PVData(final double positionSigma, final double velocitySigma) {
1377             this.positionSigma = positionSigma;
1378             this.velocitySigma = velocitySigma;
1379         }
1380 
1381     }
1382 
1383     /** Measurements types. */
1384     private static abstract class MeasurementsParser<T extends ObservedMeasurement<T>> {
1385 
1386         /** Parse the fields of a measurements line.
1387          * @param fields measurements line fields
1388          * @param stations name to stations data map
1389          * @param pvData PV measurements data
1390          * @param satRangeBias range bias due to transponder delay
1391          * @param weight base weights for measurements
1392          * @param line complete line
1393          * @param lineNumber line number
1394          * @param fileName file name
1395          * @return parsed measurement
1396          * @exception OrekitException if the fields do not represent a valid measurements line
1397          */
1398         public abstract T parseFields(String[] fields,
1399                                       Map<String, StationData> stations,
1400                                       PVData pvData,
1401                                       Bias<Range> satRangeBias, Weights weight,
1402                                       String line, int lineNumber, String fileName)
1403             throws OrekitException;
1404 
1405         /** Check the number of fields.
1406          * @param expected expected number of fields
1407          * @param fields measurements line fields
1408          * @param line complete line
1409          * @param lineNumber line number
1410          * @param fileName file name
1411          * @exception OrekitException if the number of fields does not match the expected number
1412          */
1413         protected void checkFields(final int expected, final String[] fields,
1414                                    final String line, final int lineNumber, final String fileName)
1415             throws OrekitException {
1416             if (fields.length != expected) {
1417                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
1418                                           lineNumber, fileName, line);
1419             }
1420         }
1421 
1422         /** Get the date for the line.
1423          * @param date date field
1424          * @param line complete line
1425          * @param lineNumber line number
1426          * @param fileName file name
1427          * @return parsed measurement
1428          * @exception OrekitException if the date cannot be parsed
1429          */
1430         protected AbsoluteDate getDate(final String date,
1431                                        final String line, final int lineNumber, final String fileName)
1432             throws OrekitException {
1433             try {
1434                 return new AbsoluteDate(date, TimeScalesFactory.getUTC());
1435             } catch (OrekitException oe) {
1436                 throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1437                                           "wrong date " + date +
1438                                           " at line " + lineNumber +
1439                                           " in file " + fileName +
1440                                           "\n" + line);
1441             }
1442         }
1443 
1444         /** Get the station data for the line.
1445          * @param stationName name of the station
1446          * @param stations name to stations data map
1447          * @param line complete line
1448          * @param lineNumber line number
1449          * @param fileName file name
1450          * @return parsed measurement
1451          * @exception OrekitException if the station is not known
1452          */
1453         protected StationData getStationData(final String stationName,
1454                                              final Map<String, StationData> stations,
1455                                              final String line, final int lineNumber, final String fileName)
1456             throws OrekitException {
1457             final StationData stationData = stations.get(stationName);
1458             if (stationData == null) {
1459                 throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
1460                                           "unknown station " + stationName +
1461                                           " at line " + lineNumber +
1462                                           " in file " + fileName +
1463                                           "\n" + line);
1464             }
1465             return stationData;
1466         }
1467     }
1468 
1469     /** Parser for range measurements. */
1470     private static class RangeParser extends MeasurementsParser<Range> {
1471         /** {@inheritDoc} */
1472         @Override
1473         public Range parseFields(final String[] fields,
1474                                  final Map<String, StationData> stations,
1475                                  final PVData pvData,
1476                                  final Bias<Range> satRangeBias,
1477                                  final Weights weights,
1478                                  final String line,
1479                                  final int lineNumber,
1480                                  final String fileName)
1481                                                  throws OrekitException {
1482             checkFields(4, fields, line, lineNumber, fileName);
1483             final StationData stationData = getStationData(fields[2], stations, line, lineNumber, fileName);
1484             final Range range = new Range(stationData.station,
1485                                           getDate(fields[0], line, lineNumber, fileName),
1486                                           Double.parseDouble(fields[3]) * 1000.0,
1487                                           stationData.rangeSigma,
1488                                           weights.rangeBaseWeight);
1489             if (stationData.rangeBias != null) {
1490                 range.addModifier(stationData.rangeBias);
1491             }
1492             if (satRangeBias != null) {
1493                 range.addModifier(satRangeBias);
1494             }
1495             if (stationData.rangeTroposphericCorrection != null) {
1496                 range.addModifier(stationData.rangeTroposphericCorrection);
1497             }
1498             return range;
1499         }
1500     }
1501 
1502     /** Parser for range rate measurements. */
1503     private static class RangeRateParser extends MeasurementsParser<RangeRate> {
1504         /** {@inheritDoc} */
1505         @Override
1506         public RangeRate parseFields(final String[] fields,
1507                                      final Map<String, StationData> stations,
1508                                      final PVData pvData,
1509                                      final Bias<Range> satRangeBias,
1510                                      final Weights weights,
1511                                      final String line,
1512                                      final int lineNumber,
1513                                      final String fileName)
1514                                                      throws OrekitException {
1515             checkFields(4, fields, line, lineNumber, fileName);
1516             final StationData stationData = getStationData(fields[2], stations, line, lineNumber, fileName);
1517             final RangeRate rangeRate = new RangeRate(stationData.station,
1518                                                       getDate(fields[0], line, lineNumber, fileName),
1519                                                       Double.parseDouble(fields[3]) * 1000.0,
1520                                                       stationData.rangeRateSigma,
1521                                                       weights.rangeRateBaseWeight,
1522                                                       true);
1523             if (stationData.rangeRateBias != null) {
1524                 rangeRate.addModifier(stationData.rangeRateBias);
1525             }
1526             return rangeRate;
1527         }
1528     };
1529 
1530     /** Parser for azimuth-elevation measurements. */
1531     private static class AzElParser extends MeasurementsParser<Angular> {
1532         /** {@inheritDoc} */
1533         @Override
1534         public Angular parseFields(final String[] fields,
1535                                    final Map<String, StationData> stations,
1536                                    final PVData pvData,
1537                                    final Bias<Range> satRangeBias,
1538                                    final Weights weights,
1539                                    final String line,
1540                                    final int lineNumber,
1541                                    final String fileName)
1542                                                    throws OrekitException {
1543             checkFields(5, fields, line, lineNumber, fileName);
1544             final StationData stationData = getStationData(fields[2], stations, line, lineNumber, fileName);
1545             final Angular azEl = new Angular(stationData.station,
1546                                              getDate(fields[0], line, lineNumber, fileName),
1547                                              new double[] {
1548                                                            FastMath.toRadians(Double.parseDouble(fields[3])),
1549                                                            FastMath.toRadians(Double.parseDouble(fields[4]))
1550             },
1551                                              stationData.azElSigma,
1552                                              weights.azElBaseWeight);
1553             if (stationData.refractionCorrection != null) {
1554                 azEl.addModifier(stationData.refractionCorrection);
1555             }
1556             if (stationData.azELBias != null) {
1557                 azEl.addModifier(stationData.azELBias);
1558             }
1559             return azEl;
1560         }
1561     };
1562 
1563     /** Parser for PV measurements. */
1564     private static class PVParser extends MeasurementsParser<PV> {
1565         /** {@inheritDoc} */
1566         @Override
1567         public PV parseFields(final String[] fields,
1568                               final Map<String, StationData> stations,
1569                               final PVData pvData,
1570                               final Bias<Range> satRangeBias,
1571                               final Weights weights,
1572                               final String line,
1573                               final int lineNumber,
1574                               final String fileName)
1575                                               throws OrekitException {
1576             // field 2, which corresponds to stations in other measurements, is ignored
1577             // this allows the measurements files to be columns aligned
1578             // by inserting something like "----" instead of a station name
1579             checkFields(9, fields, line, lineNumber, fileName);
1580             return new org.orekit.estimation.measurements.PV(getDate(fields[0], line, lineNumber, fileName),
1581                                                              new Vector3D(Double.parseDouble(fields[3]) * 1000.0,
1582                                                                           Double.parseDouble(fields[4]) * 1000.0,
1583                                                                           Double.parseDouble(fields[5]) * 1000.0),
1584                                                              new Vector3D(Double.parseDouble(fields[6]) * 1000.0,
1585                                                                           Double.parseDouble(fields[7]) * 1000.0,
1586                                                                           Double.parseDouble(fields[8]) * 1000.0),
1587                                                              pvData.positionSigma,
1588                                                              pvData.velocitySigma,
1589                                                              weights.pvBaseWeight);
1590         }
1591     };
1592 
1593     /** Local class for measurement-specific log.
1594      * @param T type of mesurement
1595      */
1596     private static abstract class MeasurementLog<T extends ObservedMeasurement<T>> {
1597 
1598         /** Residuals. */
1599         private final SortedSet<EstimatedMeasurement<T>> evaluations;
1600 
1601         /** Measurements name. */
1602         private final String name;
1603 
1604 
1605         /** Simple constructor.
1606          * @param name measurement name
1607          * @exception IOException if output file cannot be created
1608          */
1609         MeasurementLog(final String name) throws IOException {
1610             this.evaluations = new TreeSet<EstimatedMeasurement<T>>(new ChronologicalComparator());
1611             this.name        = name;
1612         }
1613 
1614         /** Get the measurement name.
1615          * @return measurement name
1616          */
1617         public String getName() {
1618             return name;
1619         }
1620 
1621         /** Compute residual value.
1622          * @param evaluation evaluation to consider
1623          */
1624         abstract double residual(final EstimatedMeasurement<T> evaluation);
1625 
1626         /** Add an evaluation.
1627          * @param evaluation evaluation to add
1628          */
1629         void add(final EstimatedMeasurement<T> evaluation) {
1630             evaluations.add(evaluation);
1631         }
1632 
1633         /**Create a  statistics summary
1634          */
1635         public StreamingStatistics createStatisticsSummary() {
1636             if (!evaluations.isEmpty()) {
1637                 // compute statistics
1638                 final StreamingStatistics stats = new StreamingStatistics();
1639                 for (final EstimatedMeasurement<T> evaluation : evaluations) {
1640                     stats.addValue(residual(evaluation));
1641                 }
1642                return stats;
1643 
1644             }
1645             return null;
1646         }
1647     }
1648 
1649     /** Logger for range measurements. */
1650     class RangeLog extends MeasurementLog<Range> {
1651 
1652         /** Simple constructor.
1653          * @exception IOException if output file cannot be created
1654          */
1655         RangeLog() throws IOException {
1656             super("range");
1657         }
1658 
1659         /** {@inheritDoc} */
1660         @Override
1661         double residual(final EstimatedMeasurement<Range> evaluation) {
1662             return evaluation.getEstimatedValue()[0] - evaluation.getObservedMeasurement().getObservedValue()[0];
1663         }
1664 
1665     }
1666 
1667     /** Logger for range rate measurements. */
1668     class RangeRateLog extends MeasurementLog<RangeRate> {
1669 
1670         /** Simple constructor.
1671          * @exception IOException if output file cannot be created
1672          */
1673         RangeRateLog() throws IOException {
1674             super("range-rate");
1675         }
1676 
1677         /** {@inheritDoc} */
1678         @Override
1679         double residual(final EstimatedMeasurement<RangeRate> evaluation) {
1680             return evaluation.getEstimatedValue()[0] - evaluation.getObservedMeasurement().getObservedValue()[0];
1681         }
1682 
1683     }
1684 
1685     /** Logger for azimuth measurements. */
1686     class AzimuthLog extends MeasurementLog<Angular> {
1687 
1688         /** Simple constructor.
1689          * @exception IOException if output file cannot be created
1690          */
1691         AzimuthLog() throws IOException {
1692             super("azimuth");
1693         }
1694 
1695         /** {@inheritDoc} */
1696         @Override
1697         double residual(final EstimatedMeasurement<Angular> evaluation) {
1698             return FastMath.toDegrees(evaluation.getEstimatedValue()[0] - evaluation.getObservedMeasurement().getObservedValue()[0]);
1699         }
1700 
1701     }
1702 
1703     /** Logger for elevation measurements. */
1704     class ElevationLog extends MeasurementLog<Angular> {
1705 
1706         /** Simple constructor.
1707          * @param home home directory
1708          * @param baseName output file base name (may be null)
1709          * @exception IOException if output file cannot be created
1710          */
1711         ElevationLog() throws IOException {
1712             super("elevation");
1713         }
1714 
1715         /** {@inheritDoc} */
1716         @Override
1717         double residual(final EstimatedMeasurement<Angular> evaluation) {
1718             return FastMath.toDegrees(evaluation.getEstimatedValue()[1] - evaluation.getObservedMeasurement().getObservedValue()[1]);
1719         }
1720 
1721     }
1722 
1723     /** Logger for position measurements. */
1724     class PositionLog extends MeasurementLog<PV> {
1725 
1726         /** Simple constructor.
1727          * @param home home directory
1728          * @param baseName output file base name (may be null)
1729          * @exception IOException if output file cannot be created
1730          */
1731         PositionLog() throws IOException  {
1732             super("position");
1733         }
1734 
1735         /** {@inheritDoc} */
1736         @Override
1737         double residual(final EstimatedMeasurement<PV> evaluation) {
1738             final double[] theoretical = evaluation.getEstimatedValue();
1739             final double[] observed    = evaluation.getObservedMeasurement().getObservedValue();
1740             return Vector3D.distance(new Vector3D(theoretical[0], theoretical[1], theoretical[2]),
1741                                      new Vector3D(observed[0],    observed[1],    observed[2]));
1742         }
1743 
1744     }
1745 
1746     /** Logger for velocity measurements. */
1747     class VelocityLog extends MeasurementLog<PV> {
1748 
1749         /** Simple constructor.
1750          * @param home home directory
1751          * @param baseName output file base name (may be null)
1752          * @exception IOException if output file cannot be created
1753          */
1754         VelocityLog() throws IOException {
1755             super("velocity");
1756         }
1757 
1758         /** {@inheritDoc} */
1759         @Override
1760         double residual(final EstimatedMeasurement<PV> evaluation) {
1761             final double[] theoretical = evaluation.getEstimatedValue();
1762             final double[] observed    = evaluation.getObservedMeasurement().getObservedValue();
1763             return Vector3D.distance(new Vector3D(theoretical[3], theoretical[4], theoretical[5]),
1764                                      new Vector3D(observed[3],    observed[4],    observed[5]));
1765         }
1766 
1767     }
1768 
1769     /** Input parameter keys. */
1770     private static enum ParameterKey {
1771         ORBIT_DATE,
1772         ORBIT_CIRCULAR_A,
1773         ORBIT_CIRCULAR_EX,
1774         ORBIT_CIRCULAR_EY,
1775         ORBIT_CIRCULAR_I,
1776         ORBIT_CIRCULAR_RAAN,
1777         ORBIT_CIRCULAR_ALPHA,
1778         ORBIT_EQUINOCTIAL_A,
1779         ORBIT_EQUINOCTIAL_EX,
1780         ORBIT_EQUINOCTIAL_EY,
1781         ORBIT_EQUINOCTIAL_HX,
1782         ORBIT_EQUINOCTIAL_HY,
1783         ORBIT_EQUINOCTIAL_LAMBDA,
1784         ORBIT_KEPLERIAN_A,
1785         ORBIT_KEPLERIAN_E,
1786         ORBIT_KEPLERIAN_I,
1787         ORBIT_KEPLERIAN_PA,
1788         ORBIT_KEPLERIAN_RAAN,
1789         ORBIT_KEPLERIAN_ANOMALY,
1790         ORBIT_ANGLE_TYPE,
1791         ORBIT_TLE_LINE_1,
1792         ORBIT_TLE_LINE_2,
1793         ORBIT_CARTESIAN_PX,
1794         ORBIT_CARTESIAN_PY,
1795         ORBIT_CARTESIAN_PZ,
1796         ORBIT_CARTESIAN_VX,
1797         ORBIT_CARTESIAN_VY,
1798         ORBIT_CARTESIAN_VZ,
1799         MASS,
1800         IERS_CONVENTIONS,
1801         INERTIAL_FRAME,
1802         PROPAGATOR_MIN_STEP,
1803         PROPAGATOR_MAX_STEP,
1804         PROPAGATOR_POSITION_ERROR,
1805         BODY_FRAME,
1806         BODY_EQUATORIAL_RADIUS,
1807         BODY_INVERSE_FLATTENING,
1808         CENTRAL_BODY_DEGREE,
1809         CENTRAL_BODY_ORDER,
1810         OCEAN_TIDES_DEGREE,
1811         OCEAN_TIDES_ORDER,
1812         SOLID_TIDES_SUN,
1813         SOLID_TIDES_MOON,
1814         THIRD_BODY_SUN,
1815         THIRD_BODY_MOON,
1816         DRAG,
1817         DRAG_CD,
1818         DRAG_CD_ESTIMATED,
1819         DRAG_AREA,
1820         SOLAR_RADIATION_PRESSURE,
1821         SOLAR_RADIATION_PRESSURE_CR,
1822         SOLAR_RADIATION_PRESSURE_CR_ESTIMATED,
1823         SOLAR_RADIATION_PRESSURE_AREA,
1824         GENERAL_RELATIVITY,
1825         TRANSPONDER_DELAY_BIAS,
1826         TRANSPONDER_DELAY_BIAS_MIN,
1827         TRANSPONDER_DELAY_BIAS_MAX,
1828         TRANSPONDER_DELAY_BIAS_ESTIMATED,
1829         GROUND_STATION_NAME,
1830         GROUND_STATION_LATITUDE,
1831         GROUND_STATION_LONGITUDE,
1832         GROUND_STATION_ALTITUDE,
1833         GROUND_STATION_POSITION_ESTIMATED,
1834         GROUND_STATION_RANGE_SIGMA,
1835         GROUND_STATION_RANGE_BIAS,
1836         GROUND_STATION_RANGE_BIAS_MIN,
1837         GROUND_STATION_RANGE_BIAS_MAX,
1838         GROUND_STATION_RANGE_BIAS_ESTIMATED,
1839         GROUND_STATION_RANGE_RATE_SIGMA,
1840         GROUND_STATION_RANGE_RATE_BIAS,
1841         GROUND_STATION_RANGE_RATE_BIAS_MIN,
1842         GROUND_STATION_RANGE_RATE_BIAS_MAX,
1843         GROUND_STATION_RANGE_RATE_BIAS_ESTIMATED,
1844         GROUND_STATION_AZIMUTH_SIGMA,
1845         GROUND_STATION_AZIMUTH_BIAS,
1846         GROUND_STATION_AZIMUTH_BIAS_MIN,
1847         GROUND_STATION_AZIMUTH_BIAS_MAX,
1848         GROUND_STATION_ELEVATION_SIGMA,
1849         GROUND_STATION_ELEVATION_BIAS,
1850         GROUND_STATION_ELEVATION_BIAS_MIN,
1851         GROUND_STATION_ELEVATION_BIAS_MAX,
1852         GROUND_STATION_AZ_EL_BIASES_ESTIMATED,
1853         GROUND_STATION_ELEVATION_REFRACTION_CORRECTION,
1854         GROUND_STATION_RANGE_TROPOSPHERIC_CORRECTION,
1855         GROUND_STATION_IONOSPHERIC_CORRECTION,
1856         RANGE_MEASUREMENTS_BASE_WEIGHT,
1857         RANGE_RATE_MEASUREMENTS_BASE_WEIGHT,
1858         AZIMUTH_MEASUREMENTS_BASE_WEIGHT,
1859         ELEVATION_MEASUREMENTS_BASE_WEIGHT,
1860         PV_MEASUREMENTS_BASE_WEIGHT,
1861         PV_MEASUREMENTS_POSITION_SIGMA,
1862         PV_MEASUREMENTS_VELOCITY_SIGMA,
1863         RANGE_OUTLIER_REJECTION_MULTIPLIER,
1864         RANGE_OUTLIER_REJECTION_STARTING_ITERATION,
1865         RANGE_RATE_OUTLIER_REJECTION_MULTIPLIER,
1866         RANGE_RATE_OUTLIER_REJECTION_STARTING_ITERATION,
1867         AZ_EL_OUTLIER_REJECTION_MULTIPLIER,
1868         AZ_EL_OUTLIER_REJECTION_STARTING_ITERATION,
1869         PV_OUTLIER_REJECTION_MULTIPLIER,
1870         PV_OUTLIER_REJECTION_STARTING_ITERATION,
1871         MEASUREMENTS_FILES,
1872         OUTPUT_BASE_NAME,
1873         ESTIMATOR_OPTIMIZATION_ENGINE,
1874         ESTIMATOR_LEVENBERG_MARQUARDT_INITIAL_STEP_BOUND_FACTOR,
1875         ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE,
1876         ESTIMATOR_NORMALIZED_PARAMETERS_CONVERGENCE_THRESHOLD,
1877         ESTIMATOR_MAX_ITERATIONS,
1878         ESTIMATOR_MAX_EVALUATIONS;
1879     }
1880 
1881 }